#include <TMath.h>
#include <TH1.h>
#include <TH2.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TTree.h>
#include <TBranch.h>
#include <TLorentzVector.h>

#include <iostream>
Double_t MyCosine(Double_t* mom, Double_t* ver1, Double_t* ver2);
Double_t MyDCA(Double_t* mom2, Double_t* ver1, Double_t* ver2);
Double_t MyCosine2(Double_t* mom1, Double_t* mom2);

void MyAnalyzer() {

//  ******************************************************
//             Decalaration and construction
//  ******************************************************

//  **********Variables and arrays for branches***********

   Int_t RunNumber;
   Int_t EventNumber;
   Double_t PVERTEXx;
   Double_t PVERTEXy;
   Double_t PVERTEXz;

// Lambdas
   Int_t NPAIR_L;
   Double_t DCA_L[1500];
   Double_t MASS_L[1500];
   Double_t CPA_L[1500];
   Double_t IP1_L[1500];
   Double_t IP2_L[1500];
   Double_t VERTEXx_L[1500];
   Double_t VERTEXy_L[1500];
   Double_t VERTEXz_L[1500];
   Double_t DAUGHTER1_MOMx_L[1500]; // always pos first
   Double_t DAUGHTER1_MOMy_L[1500];
   Double_t DAUGHTER1_MOMz_L[1500];
   Double_t DAUGHTER2_MOMx_L[1500];
   Double_t DAUGHTER2_MOMy_L[1500];
   Double_t DAUGHTER2_MOMz_L[1500];
   Double_t MOMENTUMx_L[1500];
   Double_t MOMENTUMy_L[1500];
   Double_t MOMENTUMz_L[1500];
   Int_t INDEX1_L[1500];
   Int_t INDEX2_L[1500];
// D mesons
   Int_t NPAIR_D;
   Double_t DCA_D[1500];
   Double_t MASS_D[1500];
   Double_t CPA_D[1500];
   Double_t IP1_D[1500];
   Double_t IP2_D[1500];
   Double_t VERTEXx_D[1500];
   Double_t VERTEXy_D[1500];
   Double_t VERTEXz_D[1500];
   Double_t DAUGHTER1_MOMx_D[1500]; // always pos first
   Double_t DAUGHTER1_MOMy_D[1500];
   Double_t DAUGHTER1_MOMz_D[1500];
   Double_t DAUGHTER2_MOMx_D[1500];
   Double_t DAUGHTER2_MOMy_D[1500];
   Double_t DAUGHTER2_MOMz_D[1500];
   Double_t MOMENTUMx_D[1500];
   Double_t MOMENTUMy_D[1500];
   Double_t MOMENTUMz_D[1500];
   Int_t INDEX1_D[1500];
   Int_t INDEX2_D[1500];

//  ****************Branches******************

   TFile *ff = new TFile("output.root");
   TTree *fTree = (TTree*)ff->Get("t1");   
   fTree->SetBranchAddress("fRunNumber", &RunNumber);
   fTree->SetBranchAddress("fEventNumber", &EventNumber);
   fTree->SetBranchAddress("fPVERTEXx", PVERTEXx);
   fTree->SetBranchAddress("fPVERTEXy", PVERTEXy);
   fTree->SetBranchAddress("fPVERTEXz", PVERTEXz);
// Lambdas
   fTree->SetBranchAddress("fNPAIR_L", &NPAIR_L);
   fTree->SetBranchAddress("fDCA_L", DCA_L);
   fTree->SetBranchAddress("fIP1_L", IP1_L);          // IP of pos, neg trki
   fTree->SetBranchAddress("fIP2_L", IP2_L);          // IP of pos, neg trk
   fTree->SetBranchAddress("fMASS_L", MASS_L);        // Inv Mass
   fTree->SetBranchAddress("fVERTEXx_L", VERTEXx_L);      // Vertex point of Ppi pair
   fTree->SetBranchAddress("fVERTEXy_L", VERTEXy_L);      // Vertex point of Ppi pair
   fTree->SetBranchAddress("fVERTEXz_L", VERTEXz_L);      // Vertex point of Ppi pair
   fTree->SetBranchAddress("fCPA_L", CPA_L);   // CPA from Primary vertex 
   fTree->SetBranchAddress("fINDEX1_L", INDEX1_L);       // Track ID of pos, neg trk :pos first
   fTree->SetBranchAddress("fINDEX2_L", INDEX2_L);       // Track ID of pos, neg trk :pos first
   fTree->SetBranchAddress("fMOMENTUMx_L", MOMENTUMx_L);    // Momentum of pair
   fTree->SetBranchAddress("fMOMENTUMy_L", MOMENTUMy_L);    // Momentum of pair
   fTree->SetBranchAddress("fMOMENTUMz_L", MOMENTUMz_L);    // Momentum of pair
   fTree->SetBranchAddress("fDAUGHTER1_MOMx_L", DAUGHTER1_MOMx_L); //Momentum of daughter tracks : Pos. trk first
   fTree->SetBranchAddress("fDAUGHTER1_MOMy_L", DAUGHTER1_MOMy_L); //Momentum of daughter tracks : Pos. trk first
   fTree->SetBranchAddress("fDAUGHTER1_MOMz_L", DAUGHTER1_MOMz_L); //Momentum of daughter tracks : Pos. trk first
   fTree->SetBranchAddress("fDAUGHTER2_MOMx_L", DAUGHTER2_MOMx_L); //Momentum of daughter tracks : Pos. trk first
   fTree->SetBranchAddress("fDAUGHTER2_MOMy_L", DAUGHTER2_MOMy_L); //Momentum of daughter tracks : Pos. trk first
   fTree->SetBranchAddress("fDAUGHTER2_MOMz_L", DAUGHTER2_MOMz_L); //Momentum of daughter tracks : Pos. trk first
// D mesons
   fTree->SetBranchAddress("fNPAIR_D", &NPAIR_D);
   fTree->SetBranchAddress("fDCA_D", DCA_D);
   fTree->SetBranchAddress("fIP1_D", IP1_D);          
   fTree->SetBranchAddress("fIP2_D", IP2_D);          
   fTree->SetBranchAddress("fMASS_D", MASS_D);        
   fTree->SetBranchAddress("fVERTEXx_D", VERTEXx_D);     
   fTree->SetBranchAddress("fVERTEXy_D", VERTEXy_D);     
   fTree->SetBranchAddress("fVERTEXz_D", VERTEXz_D);     
   fTree->SetBranchAddress("fCPA_D", CPA_D);   
   fTree->SetBranchAddress("fINDEX1_D", INDEX1_D);       
   fTree->SetBranchAddress("fINDEX2_D", INDEX2_D);       
   fTree->SetBranchAddress("fMOMENTUMx_D", MOMENTUMx_D);    
   fTree->SetBranchAddress("fMOMENTUMy_D", MOMENTUMy_D);    
   fTree->SetBranchAddress("fMOMENTUMz_D", MOMENTUMz_D);    
   fTree->SetBranchAddress("fDAUGHTER1_MOMx_D", DAUGHTER1_MOMx_D); 
   fTree->SetBranchAddress("fDAUGHTER1_MOMy_D", DAUGHTER1_MOMy_D); 
   fTree->SetBranchAddress("fDAUGHTER1_MOMz_D", DAUGHTER1_MOMz_D); 
   fTree->SetBranchAddress("fDAUGHTER2_MOMx_D", DAUGHTER2_MOMx_D); 
   fTree->SetBranchAddress("fDAUGHTER2_MOMy_D", DAUGHTER2_MOMy_D); 
   fTree->SetBranchAddress("fDAUGHTER2_MOMz_D", DAUGHTER2_MOMz_D); 

//  ************Analysis variables**************

   Int_t nentry = fTree->GetEntries();
// Int_t nentry = 1000000;
   Double_t R2_L, R2_D, R_LH;
   Double_t mom[1000][3] = {};
   Double_t ver[1000][3] = {};
   Double_t mom1[3], mom2[3], mom3[3], ver1[3], ver2[3], temp[3], pver[3];
   Double_t mass1, mass2, aaa, mom_H, lifetime;
  
   Double_t TrkID[1000][2] = {};
   Double_t mass[1000];



//  ******************************************************
//                 Analysis cuts & Constants
//  ******************************************************


	const Double_t mass_sigma = 0.0015;
   const Double_t lambdamass = 1.11568;    //defualt 1.11568
   const Double_t pionmass = TDatabasePDG::Instance()->GetParticle(kPiMinus )->Mass();
   const Double_t massrange = 4*mass_sigma;   // sigma : 0.0015  : sensitive to this
   const Double_t protonmass = TDatabasePDG::Instance()->GetParticle(kProton)->Mass();
	const Double_t kaonmass = TDatabasePDG::Instance()->GetParticle(kKaonPlus)->Mass();
// TCP Ncluster > 60 applied

//  *******************Lambda cut************************
   const Double_t DCAmax_L  = 0.15;  // defualt 1.5
   const Double_t pIPmin_L  = 0.04;  // defualt 0.05
   const Double_t nIPmin_L  = 0.04;
   const Double_t CPAmin_L  = 0.998;   // defualt 0.9
 


//  *******************D meson cut***********************
   const Double_t DCAmax_D  = 0.15;  // defualt 1.5
   const Double_t pIPmin_D  = 0.04;  // defualt 0.05
   const Double_t nIPmin_D  = 0.04;



//  *******************Theta_cs cut**********************
//   const Double_t CPAmin_H = 0.98;   // default 0.98
//   const Double_t CosineMin = 0.98;  // cosine angle between lambda and p_pi pair
//   const Double_t Rmin_H = 0.5;   // 1.2 
//   const Double_t IPmin_H = 0.05;


// Out file
   TFile *fout = new TFile("result.root","RECREATE");   

//  ****************************************************	
//                  Main loop start
//  ****************************************************   

   for (Int_t i=0;i<nentry;i++) { 
      fTree->GetEntry(i);
      if (i%1000000==0) cout  << (Double_t)i/nentry*100 << "%" << " for H0"  << endl;
//  *******************Track cuts**********************
     for(Int_t j=0;j<NPAIR_L;j++) { // Lambda loop
		 if (NPAIR_L ==0) continue;  // reject 0 event
	      if (DCA[j] > DCAmax_L) continue;
	      if (TMath::Abs(IP1[j]) < pIPmin_L || TMath::Abs(IP2[j]) < nIPmin_L)continue;
       	if (TMath::Abs(MASS[j] - lambdamass) > massrange_L) continue;
// Cuts for Lambdas : CPA cut need 


     for(Int_t l=0;l<NPAIR_D;l++) { // K Pi loop
	      if (DCA[l] > DCAmax_D) continue;
	      if (TMath::Abs(IP1[l]) < pIPmin_D || TMath::Abs(IP2[l]) < nIPmin_D)continue;
       	if (TMath::Abs(MASS[l] - Dmass) > massrange_D) continue;
// **************need more cuts and resolution dist****


// under construction.....


//  ************Extracting physics value***************
       R2_L = (PVERTEXx - VERTEXx[j])*(PVERTEXx - VERTEXx[j]) + (PVERTEXy - VERTEXy[j])*(PVERTEXy - VERTEXy[j]) + (PVERTEXz - VERTEXz[j])*(PVERTEXz - VERTEXz[j]);
       ver1[0] = VERTEXx[j];  ver1[1] = VERTEXy[j];  ver1[2] = VERTEXz[j]; // I set Lambda have index 1      j : 1,   l : 2
	    pver[0] = PVERTEXx;  pver[1] = PVERTEXy;  pver[2] = PVERTEXz;
	    ver2[0] = VERTEXx[l];  ver2[1] = VERTEXy[l];  ver2[2] = VERTEXz[l]; // index j is outter one  : 2
       mom1[0] = MOMENTUMx[j];  mom1[1] = MOMENTUMy[j];  mom1[2] = MOMENTUMz[j];
       mom2[0] = MOMENTUMx[l];  mom2[1] = MOMENTUMy[l];  mom2[2] = MOMENTUMz[l];
       mom3[0] = mom1[0]+mom2[0];  mom3[1] = mom1[1]+mom2[1];  mom3[2] = mom1[2]+mom2[2];        mass2 = MASS[l];  mass1 = MASS[j];

/////////////////////////////////////// 
 
       if (TMath::Abs(mass2-Dmass) < massrange_D) continue;  // outter pair should be lambda
//       if (mass1-lambdamass > -0.0086) continue;   // estimation max ppi mass from MC
       if (mass1-lambdamass > -0.00568) continue;


//  **************Filling histograms*****************
       hcpaL->Fill(MyCosine(mom2, ver1, ver2));  //cut distribution
       hcpaH->Fill(MyCosine(mom3, pver, ver1));
       hdcaH->Fill(MyDCA(mom2, ver1, ver2));
       hdcaHvsR_H->Fill(TMath::Sqrt(r2_1), MyDCA(mom2, ver1, ver2));

       hR->Fill(TMath::Sqrt(r2_1));
       hR->Fill(TMath::Sqrt(r2_2));
       hangleH->Fill(MyCosine2(mom3, mom2));   
 
       mass2 = lambdamass; // I use correct lambda mass for mass2 to clean mass broadening from resolution
       mom_H = TMath::Sqrt(mom3[0]*mom3[0] + mom3[1]*mom3[1] + mom3[2]*mom3[2]);
       Double_t e_1 = TMath::Sqrt(mom1[0]*mom1[0] + mom1[1]*mom1[1] + mom1[2]*mom1[2] + mass1*mass1);
       Double_t e_2 = TMath::Sqrt(mom2[0]*mom2[0] + mom2[1]*mom2[1] + mom2[2]*mom2[2] + mass2*mass2);
       Double_t H_mass = TMath::Sqrt(mass1*mass1 + mass2*mass2 + 2*(e_1*e_2 - mom1[0]*mom2[0] - mom1[1]*mom2[1] - mom1[2]*mom2[2]));
       lifetime = (TMath::Sqrt(r2_1))*H_mass/mom_H;
       h5->Fill(H_mass);

       Double_t e11 = -999;
       Double_t e22 = -999; 
       
       if(TMath::Sqrt( mom3[0]*mom3[0] +  mom3[1]*mom3[1] + mom3[2]*mom3[2] < e_1 + e_2) ) {
       proton_4mom.Boost(-1*mom3[0]/(e_1 + e_2), -1*mom3[1]/(e_1 + e_2), -1*mom3[2]/(e_1 + e_2));
       pion_4mom.Boost(-1*mom3[0]/(e_1 + e_2), -1*mom3[1]/(e_1 + e_2), -1*mom3[2]/(e_1 + e_2));
       e11 = proton_4mom.E();
       e22 = pion_4mom.E();
       dalitz->Fill(e22,e11);
       }
        
        
       if (TMath::Sqrt(r2_1) > 0.5) hdcaH_1->Fill(MyDCA(mom2, ver1, ver2));
       if (TMath::Sqrt(r2_1) > 1) hdcaH_2->Fill(MyDCA(mom2, ver1, ver2));
       if (TMath::Sqrt(r2_1) > 1.5) hdcaH_3->Fill(MyDCA(mom2, ver1, ver2));



       if (MyDCA(mom2, ver1, ver2) > DCAmax_H) continue;    // dca cut PPi vertex and lambda track
       hangleH_afterDCA->Fill(MyCosine2(mom3, mom2));
       AngleHLvsM_H_0->Fill(MyCosine2(mom3, mom2) ,H_mass);

       if (MyCosine(mom2, ver1, ver2) < CPAmin_H) continue; // lambda cpa cut
       if (MyCosine(mom3, pver, ver1) < CPAmin_H) continue; // H cpa cut
       h9->Fill(H_mass);
       if (MyCosine2(mom3, mom2) < CosineMin) continue;     // momentum H, L angle cut
       if (TMath::Sqrt(r2_1) < Rmin_H) continue;
       h7->Fill(H_mass);
       dalitz_signal->Fill(e22,e11);




       } // l loop
     }  // j loop
   }  // i loop     



cout << "The number of events : " << nentry << endl;


TCanvas *cc = new TCanvas();
cc->Divide(2,2);
cc->cd(1);
h1_10->Draw();
cc->cd(2);
h2->Draw();
cc->cd(3);
h4->Draw();
cc->cd(4);
//h11->Draw();

fout->Write();
}  // main



Double_t MyCosine(Double_t* mom, Double_t* ver1, Double_t* ver2)  // for CPA
{
   Double_t vector[3] ={};
   vector[0] = ver2[0] - ver1[0];
   vector[1] = ver2[1] - ver1[1];
   vector[2] = ver2[2] - ver1[2];
   Double_t cosine = (vector[0]*mom[0] + vector[1]*mom[1] + vector[2]*mom[2]) / TMath::Sqrt( (vector[0]*vector[0] + vector[1]*vector[1] + vector[2]*vector[2])*(mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2]));
   return cosine;
}

Double_t MyDCA(Double_t* mom2, Double_t* ver1, Double_t* ver2)  // outter one is ver2, mom2
{
   Double_t t;
   Double_t p_2 = mom2[0]*mom2[0] + mom2[1]*mom2[1] + mom2[2]*mom2[2];
   Double_t pv = mom2[0]*(ver1[0]-ver2[0]) + mom2[1]*(ver1[1]-ver2[1]) + mom2[2]*(ver1[2]-ver2[2]);
   if (p_2 == 0) return -999;
   t = pv/p_2;
   Double_t r2 = (mom2[0]*t + ver2[0] - ver1[0])*(mom2[0]*t + ver2[0] - ver1[0]) + (mom2[1]*t + ver2[1] - ver1[1])*(mom2[1]*t + ver2[1] - ver1[1]) + (mom2[2]*t + ver2[2] - ver1[2])*(mom2[2]*t + ver2[2] - ver1[2]);
   return TMath::Sqrt(r2);
}
Double_t MyCosine2(Double_t* mom1, Double_t* mom2)
{
   Double_t p1 = TMath::Sqrt(mom1[0]*mom1[0] + mom1[1]*mom1[1] + mom1[2]*mom1[2]); 
   Double_t p2 = TMath::Sqrt(mom2[0]*mom2[0] + mom2[1]*mom2[1] + mom2[2]*mom2[2]);
   Double_t p1p2 = mom1[0]*mom2[0] + mom1[1]*mom2[1] + mom1[2]*mom2[2]; 
   return p1p2/p1/p2;
}


